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It is shown how to apply the Maximum Entropy Method (MEM) to numerical 

Dyson-Schwinger studies for the extraction of spectral functions of correlators from 

| their corresponding Euclidean propagators. Differences to the application in lattice 

£N) ■ QCD are emphasized and, as an example, the spectral functions of massless quarks 

^ . in cold and dense matter are presented. 
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1 Introduction 



For the investigation of QCD in the strongly coupled regime, non-perturbative 
numerical methods such as lattice QCD and truncated Dyson-Schwinger equa- 
tions can be employed. Similar to lattice QCD, numerical Dyson-Schwinger 
studies of QCD n-point functions are usually performed in Euclidean space (for 
recent reviews see [1,2,3]). Although Dyson-Schwinger studies rely on trunca- 
tion schemes, they have the great advantage, that they can in principle be 
solved in the continuum limit and with much higher numerical accuracy. 

An extraction of dynamical properties, in particular the spectral functions of 
propagators, is however highly desirable. In quantum Monte Carlo simulations, 
the Maximum Entropy Method (MEM) turns out to be an especially suited 
tool and has been successfully applied in condensed matter physics (see [4] 
for a review), lattice QCD in the vacuum (see [5] for a review), as well as at 
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finite temperatures [6]. In this work we show, that it can also successfully be 
employed for the extraction of spectral functions in Dyson- Schwinger studies. 

The starting point for the MEM is the linear relation between the spectral 
function and numerically determined Euclidean correlation functions via gen- 
eralized Kallen-Lehmann representations. Since the inversion of the latter is 
in general ill-posed due to the spectral properties of the linear operator, fur- 
ther knowledge has to be implemented non-linearly. This can be done from a 
regularization point of view leading to the "historical maximum entropy" [7] 
or in a Bayesian approach leading to the "classic maximum entropy" [8,9] and 
to "Bryan's method" [10]. 

The key idea in the latter case is the interpretation of the spectral function as a 
probability distribution due to its special properties and a proper consideration 
of the numerical error. As a result, the MEM determines the most plausible or 
expected spectral function for a given Euclidean correlator with known errors 
and some prior knowledge. It does not rely on a special form of the function 
and should, with decreasing errors, converge towards the exact solution. 

For numerical Dyson- Schwinger studies, the method seems especially reliable, 
since the calculations can usually be performed with much higher numerical 
accuracy and for much more momentum points of the correlators than in 
lattice QCD. 

This paper is organized as follows: In section 2 we collect those properties 
of spectral functions for fermions and bosons, which make them suitable for 
the application of the MEM. In section 3 we discuss the MEM procedure 
itself, in particular its adaptation to Dyson-Schwinger studies and to fermions. 
In section 4 we sketch the numerical implementation. After this we show in 
section 5 some results for massless color-superconducting fermions in dense 
quark-gluon matter. Finally we summarize and conclude in section 6. 



2 Spectral functions and their properties 

Performing calculations for propagators in Euclidean space or rather within 
the imaginary time formalism means the determination of Matsubara prop- 
agators. These are related to spectral functions in Minkowski space via the 
generalized Kallen-Lehmann representation (see e.g. [11]). Using the Euclidean 
conventions of [1] we have 




(1) 
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for a Dirac fermion propagator, where p 4 = u n denotes a given Matsubara 
frequency and p is the chemical potential. The spectral function is given by 



p(^,p^ = ^E e ^ (Ei_/liVi) (l + e-^)(^ak)H^|Ox 

Z yl J ) l,n 

x S (w + E x - E n )S 3 (p + pi- p n ) . (2) 



Therefore p(u, p)^ is hermitian and has only positive eigenvalues in a given 
Hilbert space 1 . Furthermore it has to fulfill the sum rule 



l4 = 7T- I duJ P(^,p), (3) 



2vr 



■oo 



as a consequence of the (anti-)commutation relations. Thus p(u;,p)74/4 can 
be identified with a probability distribution, which is the key property for 
motivating the use of the MEM. 

For a massive relativistic fermion in an isotropic, even parity and T-symmetric 
phase, we can parameterize 

p(uj,p)=27r(uj-f 4 p e (uj,p) -ip-^p v (uj,p) + p s (uj,p)) (4) 
and knowing that all eigenvalues have to be positive, we get 

UJpe(uJ,P) > \fp 2 p v (uj,p) 2 + p s (uJ,P) 2 > (5) 

and furthermore the sum rules 

/oo 
dw ujp e (uj,p), 
-oo 

/oo 
dw p v {uo,p), 
-oo 

POO 

0= I dw p s (uj,p). (6) 



oc 



For the application discussed in section 5, we will for simplicity restrict ourself 
to the chiral limit, i.e. massless fermions (the case of massive fermions is 
discussed at the end of section 3). In this case, we can rewrite the propagators 
by using the energy projectors A 1 * 1 = | (l ± ^74^) an d obtain 

1 This argument does therefore not hold e.g. for gauge fields. 
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p(uj,p) = 2tt (p+(w,p)A + 74 + p (w,p)A 74) 
The spectral functions p^ then fulfill 



(7) 



p ± (u,p)>0, 



/oo 
du p ± (uj,p) = 1, 
-00 



(8) 



with 



-00 — «p 4 + /! — UJ 



(9) 



For completeness we also wish to show, how the solutions of the inhomogeneous 
Bethe-Salpeter equation (BSE) can be used to determine the spectral functions 
of mesons or diquarks. For a given current J a (x) = ^(x)T a ijj(x), the solution 
for the corresponding Bethe-Salpeter amplitude T a (q;P) in momentum space 
(see e.g. [1]) determines the time-ordered product 



x s( q + ^)V(i;P)S(q-^), (10) 



such that the current-current correlator is given by 



(T J"(x)J b (y))e = Mm TY (<T J°(x)«i( !/ ),H^>^') 

with 



= (S(q + ^)T a (q; P)S(q - ^)T h ) . (12) 

Again D ab (P) possesses a generalized Kallen-Lehmann representation and the 
spectral function has to fulfill positivity conditions. 
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3 Maximum Entropy Method (MEM) 



The Maximum Entropy Method is a numerical tool for the inversion of poten- 
tially ill-posed linear equations by the implementation of additional informa- 
tion, i.e. constraints. It can be easily viewed from a standpoint of regulariza- 
tion as adding some non-linear auxiliary conditions, leading to the so-called 
"historical maximum entropy" [7]. This is known to underfit the data by over- 
estimating the effective number of degrees of freedom, thus leading to solutions 
that are closer to the prior estimate. Usually, the somewhat converse Bayesian 
viewpoint is considered, since it essentially adjusts the number of effective de- 
grees of freedom and also allows for an error estimation [4,5]. We briefly review 
this here, emphasizing the adaption to our problem. 

Given a (numerically evaluated) Euclidean correlator, which is treated as 
'data' D, the objective is to determine the most plausible (related to "classic 
maximum entropy" [8,9]) or the most expected (related to "Bryan's method" [10]) 
spectral function Pmem by taking into account prior knowledge H{m) of the 
solution, regulated by the prior estimate m to be defined below. The key 
entity is the plausibility functional P[p\DH(m)} for the spectral function p 
under given D and H(m). With help of the so called "Bayesian theorem" for 
conditional plausibilities 

P[XY] = P[X\Y]P[Y] = P[Y\X]P[X], (13) 

this can be brought into the form 



mDH(m)] = Wp?MW^ 

[Pl [ n P[D\H{m)\ 

oc P[D\pH(m)]P[p\H(m)}. (14) 

Here, we have introduced the "likelihood function" P[D\pH(m)] for the plau- 
sibility of the data D under given p and if (to) and the "prior probability" 
P[p\H(m)} for the plausibility of p under the prior knowledge if (to). The con- 
stant plausibility P[D\H(m)] of the data D under the prior knowledge H(m) 
can be dropped, since we normalize the plausibility functional at the end. 

Considering the function (or sequence at finite temperature) D in an interval 
[a, b] as uncorrected data points obeying a Gaussian distribution functional 2 , 
we have 



P[D\pH(m)]=exp(-L\p]), (15) 



The justification of this will be discussed for a given application. 



5 



with the likelihood 



L[p] 



b — a 



b — a 



1 



1 



J a 




'* 2(7? ' 



(16) 



where D[p] is given by the generalized Kallen-Lehmann representation and a 
measure VD for the discretized integral 



The prior probability P[p\H{m)\ is usually somewhat arbitrary and essentially 
implements the positivity conditions non-linearly, at least from the regular- 
ization point of view. It can be motivated by the law of large numbers or 
axiomatically constructed (see e.g. [4,5]). The key idea in the latter case is 
to consider the spectral function as a probability distribution and derive the 
most general functional, fulfilling the requirements of subset independence, 
coordinate invariance, system independence and scaling. It is then of the form 



with the scaling factor a for the prior estimate m and the plausibility for the 
scaled prior estimate 




(17) 




(18) 




(19) 



where S [p] is the negative semi-definite entropy 




(20) 



with the measure 




(21) 
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in the saddle point approximation around the prior estimate. The scaling 
factor a basically scales the maximum of the entropy and will be balanced by 
the likelihood. Furthermore we will assume, that the plausibility P[a\H{m)\ 
can be dropped, which is called Laplace's rule and can be justified a posteriori, 
as discussed in our example in section 5. 

From Eq. (14) we therefore finally get 



with the negative semi-definite functional Q [p] = aS [p] —L[p], Z determined 
by normalization and the measure 



It is worth noting, that the "historic maximum entropy" [7] simply determines 
the maximum, which is unique if it exists [5], of the functional Q [p] with a 
chosen, such that L — 1. We will however follow "Bryan's method" [10], aiming 
at the most expected spectral function, by computing 



where it is assumed that P[p\DH{m)} is sharply peaked around its maximum 
p a . We therefore define 




(22) 




(23) 




(24) 




(25) 



with {Afe} being the eigenvalues of 




(26) 



and finally get the most expected spectral function via 
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Pmem — 



Jo 



/ da p a P[a\DH(m)] 



(27) 



It should be noted, that P[a\DH(m)] is formally not integrable due to the 
saddle point approximation. However, this becomes only relevant for very 
large values of a. In any numerically considered interval the function decreases 
exponentially for precise enough and many data points. For consistency, the 
choice for the upper cutoff for a should always be quoted. In the "classical 
maximum entropy" , the most plausible spectral function p a is determined by 
maximizing Q [p] and P[a\DH{m)\ simultaneously. In our case, due to (in 
principle) arbitrarily many data points, this turns out and is known [4] to 
agree with the most expected spectral function Pmem- 

In comparison to previous applications, we have formulated the MEM for arbi- 
trarily discretized functions D and p a , since we want to deal with (in principle) 
continuous functions from truncated Dyson-Schwinger calculations. Therefore, 
the single- value decomposition as proposed in the Bryan algorithm [10] for the 
numerical determination of p a does not work. However, our new treatment of 
the spectral function opens the possibility of a better suited discretization, 
which can be adopted to the specific form of p a . In this way, we are also 
able to significantly reduce the number of points, which are needed for the 
numerically discretized spectral function. 

As already mentioned, the Bayesian approach offers the possibility of an error 
estimation. If we consider an interval I = [uj 1 ,uj2] of the spectral function, the 
expectation value of p a for fixed a in this interval for constant weighting is 
given by 




(28) 



and, hence, for the most plausible function by 




(29) 



For the variance around this value, we first consider [8,9] 



{6p a (ui)6p a (uj))~4: 





d 2 Q 



(30) 



dpidpj 



P=Pa 
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with fixed a and 5p a (u) = piu) — p a {uj). Therefore in the given interval /, we 
get 



{u) 2 — OJl) 2 Jlxl 



dudu' 



5 2 Q 



(31) 



P=Pa 



and 



((5puEu) 2 )i = J™da ((5p Q ) 2 ) j P[a|M(m)]. (32) 

At the end of this section, we want to indicate, how the MEM can be applied 
for massive fermions. For p e (u,p), the upper procedure can be performed 
analogously due to the properties given in Eqs. (5) and Eqs.(6). On the other 
hand, we would need to extend the procedure for the whole propagator by 
utilizing that p(u,p) , ~f4 1 /4 can be considered as a probability distribution. With 
the propagator again denoted as data D, we generalize 



1 fb Tr {{Dip,) - D[p] (p 4 )) t (D(p 4 ) - D[p](p 4 ))) 
L[p] ^—a L dPi W ' 

/OO 
du Tr ((p(u;) - m(uj) - In (p(cj)m(cu)" 1 ) p(u)) 74) , (33) 

where the prior estimate m is now matrix- valued. Since p(u,p)'~f4/4 is hermi- 
tian and positive, it can be written as g'ppg, with g e U(4) and po a diagonal 
matrix with positive eigenvalues. The unconstrained measure T>p then becomes 



Vp^da\[ J 2aA ^ dXl) (34) 

i,a ' ^ 

where p ijCl is the a-th eigenvalue and d\ the Haar measure for the group 
U(4) at an energy uji. However, the path integral is usually constrained by 
symmetries, i.e. p(a;,p)74 = w p(u,p)'j4h for h e Ti,. It can be easily seen that 
the group needs to be only integrated over the factor group of U(4) and the 
conjugate closure of TC and that only the independent eigenvalues need to be 
considered. For the case given by Eq. (4), we obtain with 



p(u,p)=p + (iu,p)A+^ 4 + p (u,p) K,pli » ( 35 ) 
where A J- = \{l ± (i cos 6 (uj,p)-f 4 ^ - sin e(u,p)~/ 4 )), simply 
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v P - da n ^/ 2aA ^ ft+ J 2aAa ^~ ^. 

X . X V 7T V VT 7T 



In the approximation of setting the integrand of the ^-integration to be con- 
stant and equal to the value 6i given by the maximization of the functional 
Q[p], the practical MEM procedure again becomes similar to the upper case. 



4 Numerical implementation 



For the numerical determination of the most plausible spectral function, we 
first need the data D with errors a on an interval [a, b\. Furthermore we have 
to choose a suitable interval [u;i,u;2] for the considered part of the spectral 
function and a prior estimate m. The interval [u^u^] is usually suggested by 
the involved scales and can be chosen to be rather large. The prior estimate 
m is in our case taken as a constant and can be estimated from the sum rules 
(Eqs. (6)), if the spectral function varies only in a certain interval AI. It can 
also be adopted to the knowledge from other methods. Since the main purpose 
is the implementation of positivity, the results turn out to be comparatively 
insensitive to its choice. The advantage of our method is that, after an eventual 
test calculation, we can choose the discretization, i.e. the abscissas and weights, 
of the spectral function, such that it is interpolated by a small number of 
points. 



We summarize the procedure as follows: 



• Take the data D with errors a on an interval [a, b] and choose a prior 
estimate m on an interval [c^i,^]- 

• Determine the maximum of the functional Q[p] = L[p] — aS[p] for a fixed 
value of a. Due to its simple form, the Marquardt-Levenberg method [7] is 
very well suited. 

• Choose a discretized interval for a, such that P[ct\DH(m)] is strongly 
peaked. Z is determined by normalization. For consistency, the choice for 
the upper cutoff for a should always be quoted. 

• Calculate pmem and eventually (pmem)j and ((^Pmem) 2 )/ for a chosen in- 
terval I. 



10 



5 Spectral functions of quarks in cold dense matter 

5.1 Color- superconducting quark matter 

As an example, we now want to present results for spectral functions of mass- 
less quarks in dense matter at vanishing temperature, as they have been deter- 
mined in [12]. We consider the gapped channel in the color-superconducting 
2SC phase at a quark chemical potential of /i — lGeV. The propagator is 
then of the form 



and S + (p^p) is related to p + (u,p) by Eq. (9). For the bare normal quark 
propagator with a constant gap A, the spectral function is then given by 



with p — \p\ and — \J (p — //) + A 2 . We will see below, how a non-trivial 
p 4 -dependence generates a finite width. 

5.2 Input data and error estimate 

As described in section 3, the main input for the MEM is the data with a 
proper error estimate. In Dyson-Schwinger studies, these are obtained by self- 
consistent solutions of truncated integral equations. To lowest order, the error 
of S + therefore scales with the error of the numerical integrals, which deter- 
mine the normal and anomal self energies S + and $ + (see [12] for details). In 
our case, we have chosen a simple Riemann quadrature for the multidimen- 
sional integrals, due to the principle-value-type behavior around the Fermi 
surface for ungapped channels. The error is therefore of order 0(h), where h 
is a scaling factor of the integration mesh. For the error estimation, we there- 
fore calculate the propagator for two different h and extrapolate linearly to 
h — 0. The data are then taken as the result for the smaller scaling factor h 
and the error as the difference between these data and the extrapolated result. 
In addition, the errors around nearest neighbors are averaged in order to avoid 
(artifically appearing) vanishing errors. 

We also have to justify, that correlations between the data points are negligible 



S(p4,p) = S + (p 4 ,p)A^ 4 + S {P4,P)^14 



(37) 




(38) 




11 




p 4 [GeV] 

Figure 1. Real and imaginary part of the quasiparticle propagator S + in the gapped 
channel of the 2SC phase at a chemical potential of \i = lGeV and momentum 
p = 0.9GeV as a function of the Euclidean energy The errors are given as 
shaded regions around the lines and are of the order of their thickness. 



for the likelihood in the form given in Eq. (16). Since our numerical integrals 
for different values of p^ are in principle independent, this is assumed to be 
true, at least when the discretized data is coarser than the numerical integral 
of the self energies. 

The input data for S + {p^p) for the following example is chosen on an interval 
[OGeV, lGeV] for the gapped channel in the 2SC phase at a quark chemical 
potential of fi = lGeV. As an illustration, the input for momentum p = 
0.9GeV is shown in Fig. 1. We consider the input as continuous due to our 
many data points and it has small errors of less than 1% in absolute value 
above 30MeV. 



5.3 Choice of the prior estimate 



For the data input with given errors, we now need to choose an interval for 
the spectral function and a non-vanishing prior estimate. We take the com- 
paratively large interval [— 1.5GeV, 2GeV]. Furthermore we choose an interval 
for P[a\DH(m)] as discussed in the following subsection. For different prior 
estimates m = 0.001,0.01,0.1 and 1.0 GeV -1 at momentum p = 0.9GeV, we 
then obtain the most expected spectral functions as shown in Fig. 2. It turns 
out, that the extracted spectral function is remarkably insensitive on the vari- 
ation of the prior estimate, even when varying it by more than three orders of 
magnitude. This is mainly related to the small errors of our data. We therefore 
fix m = O.lGeV - in the following. 
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<D[GeV] co[GeV] 



Figure 2. The most expected spectral function p + (u>) in the gapped channel 
of the 2SC phase at momentum p = 0.9GeV for constant prior estimates 
m = 0.001,0.01,0.1 and 1.0 GeV -1 in linear (left) and logarithmic (right) presenta- 
tion. 

5.4 The a -dependence 

It is also interesting to investigate the a-dependence of the maximum of 
the functional Q[p) as well as of P[a\DH{m)}. Since P[a\DH{m)] shows a 
pronounced maximum at a max , we choose the interval / = [ai ow , ahigh] for 
P[a\DH(m)} to be non-vanishing and normalized, such that P[a\DH(m)} > 
KT 1 x P[a max \DH{m)} for a E P 

Again, for momentum p = 0.9GeV, the results are shown in Fig. 3. On the 
left-hand side, we show P[a\DH(m)], normalized on J. On the right-hand 
side, we present the maximum of Q[p] for a = a m i n , a max and cthigh- They are 
only weakly varying, even when comparing the border of the interval / to the 
maximum. This also substantiates Laplace's rule for P[a\H(m)] a posteriori, 
assuming that its a-dependence is weaker. 

Apart from this, in Fig. 3, the most expected spectral function is equal to the 
most probable spectral function given at a max and therefore not shown. Thus 
the "classic maximum entropy" gives very similar results as "Bryan's method" 
in our case of "many" data points and small errors. 

5. 5 Error estimate and final result 

Finally we are able to perform an error estimation around the expected value 
of a given interval. Since we have two pronounced peaks for quasiparticles and 
quasiparticle-holes in the spectral function, we choose the intervals associated 
to their full width at half maximum (FWHM). The result for momentum 
p = 0.9GeV is shown on the right in Fig. 4. Again, the errors turn out to be 
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0.9 1 1.1 

co[GeV] 



Fi gure 3. The function P[ct\DH(rn^ normalized between ct m in and otjiigh 

(see text) 

with maximum a max (left) and the maxima Pajow /max/high °f the functional Q\p] 
for given a (right). Both for the gapped channel of the 2SC phase at momentum 
p = 0.9GeV. 




Figure 4. A contour plot of most expected spectral density of the gapped 2SC phase 
at [i = lGeV as described in the text (left) and the spectral function for momentum 
p = 0.9GeV with the expectation value within the FWHM and its error estimate 
as shaded background (right). 



very small. 



On the left-hand side we show a contour plot of the spectral density as a 
function of the energy u and momentum p. The light (online yellow) line shows 
the maxima of the quasiparticle and quasiparticle-hole branches. For fixed 
momentum p, the difference between the dark (online blue) lines below and 
above the light (online yellow) line gives the FWHM of the corresponding peak. 
The latter is neglected in BCS-type spectral functions as given in Eq. (38). 
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6 Summary and conclusions 



In this paper we have outlined, how the MEM can be adapted to numerical 
Dyson-Schwinger studies, in particular to fermions. It turns out, that the ex- 
tracted spectral functions are much more reliable and stable against variation, 
than in applications of the MEM in lattice QCD. Reasons for this are the com- 
paratively small errors on the input functions and the almost arbitrary large 
number of data points. Comparing however to the systematic error, coming 
from necessary truncations in Dyson-Schwinger studies, this error is negligi- 
ble. Therefore this method can be useful for further applications in mesons or 
diquarks investigations, since currently all calculations have to be extended 
to complex momenta (see [13,14]). Avoiding this and reducing the numerical 
effort drastically, the MEM might help to improve or extent known truncation 
schemes. 
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